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Self-overlap as a method of analysis in Ising models 
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The damage spreading method (DS) provided a useful tool to obtain analytical results of the 
thermodynamics and stability of the 2D Ising model -amongst many others-, but it suffered both 
from ambiguities in its results and from large computational costs. In this paper we propose an 
alternative method, the so called self-overlap method, based on the study of correlation functions 
measured at subsequent time steps as the system evolves towards its equilibrium. Applying marko- 
vian and mean field approximations to a 2D Ising system we obtain both analytical and numerical 
results on the thermodynamics that agree with the expected behavior. We also provide some ana- 
lytical results on the stability of the system. Since only a single replica of the system needs to be 
studied, this method would seem to be free from the ambiguities that affiicted DS. It also seems to 
be numerically more efficient and analytically simpler. 

PACS numbers: 05.10-a, 05.20-y, 64.60-Ht 



I. INTRODUCTION 

The damage spreading (DS) method 1] is a remarkable 
tool amongst the many ones developed in recent years in 
the effort to understand the dynamics of cooperative sys- 
tems. Very roughly speaking the goal of the method is 
to study the stability of a cooperative system under a 
small perturbation: if perturbations die off after some 
time then the system must be in a stable, ordered state; 
if small perturbations always get amplified however then 
the system must be in a disordered, chaotic state. By 
studying how far away the final states are from the initial 
ones given the initial perturbation one can get informa- 
tion about the system such as, for example, the Lyapunov 
exponents. 

Of course there is a key aspect that differentiates 
cooperative systems from classical dynamical systems. 
Namely, that in the former case given the complexity of 
the systems under study we almost never have at our 
disposal a detailed analytical solution to the equations 
of motion in order to study the system's stability under 
perturbations. Here is where the DS method comes in: 
its operational side amounts to an algorithm designed 
to study how small perturbations spread within the 
system by working in detail how each of the system's 
components react to the changes. The method has been 
applied to many different dynamical systems such as 
Ising systems [2, 3, 4, 5, jj, 9, 10, 11], Kauffman networks 
[H, 111, El, spin glasses [ll, ll6[, cellular automata [13 
amongst others, yielding in many cases useful infor- 
mation about their evolution and stability. Succinctly 
speaking, the algorithm analyzes the evolution of two 
almost identical states of the system. The damage 
(difference between the two initial states) is specified as 
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part of the initial conditions. That is, on one side we 
have a specified state of the system, and on the other we 
have a replica that only differs in a small perturbation 
(the damage) from this original state. One then fixes 
the stochastic evolution to be the same for each replica 
(in a Monte Carlo simulation the method imposes the 
same random numbers at each step of time on both 
copies for instance). As we let the two copies evolve, the 
method analyzes their distance (Hamming distance) as 
a function of time. Useful information about the system 
can then be extracted from this, not only numerically 
but in some cases also analytically. 

However, as was shown in 0, i and (and refer- 
ences therein) DS has been shown to be ill-defined in 
the sense that different -and equally legitimate- algo- 
rithmic implementations of the same physical system's 
dynamics can yield different DS properties. This am- 
biguity stems from the fact that while the transfer ma- 
trix for the evolution of a single system is completely de- 
termined by the one-point correlation functions [9|, the 
simultaneous evolution of two replicas however is gov- 
erned by a joint transfer matrix determined by two-point 
correlation functions. For example, Glauber and both 
standard and uncorrelated heat bath (HB) algorithms 
satisfy detailed balance with respect to the same Hamil- 
tonian. It follows that these three different update rules 
generate the same equilibrium ensemble and are there- 
fore equally legitimate to mimic the evolution in time of 
an Ising system coupled to a thermal reservoir. Accord- 
ingly, the one-point correlation functions for the three 
cases coincide and the corresponding transfer matrices 
for single systems are identical. On the other hand the 
two-point functions for HB and Glauber dynamics are 
different; hence damage evolves differently in either case 
(see ,9] and references for a extended quantitative version 
of this argument). As long as the results depend on the 
algorithm being implemented, one can not assert that the 
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results obtained from a given DS analysis are conclusive 
and unambiguous. 

This handicap is a major motivation in order to 
search for an alternative method of stability analysis. 
Our goal in this paper will be to propose a different 
approach to study the stability of cooperative sys- 
tems. By relying heavily on the above mentioned fact 
that the evolution of a single system is determined 
only by the one-point correlation functions we will try 
to eliminate some of ambiguities found in the DS method. 

In order to be specific, as a test case we will focus on 
the study of a well-known type of system: Ising models. 
In 0] Votja tackled the 2D Glauber-Ising model via DS. 
He obtained results on the thermodynamics (magneti- 
zation, ferromagnetic transition) and stability (regular 
vs. chaotic behavior) of the model both analytically 
and numerically. Due to the nature of the method 
however (at every step of time we must keep account of 
the two replicas), there is obvious room for improving 
the computational efficiency. This is also the case in 
the analytical realm, where accounting for the way in 
which at each step of time the differences between the 
two copies may increase inevitably leads to lengthy 
computations (as was shown in 0, H, 0, 0, @1)- This on 
itself constitutes a second motivation in order to search 
for alternative methods. 



The method of analysis that we will prop ose here, 
so-called the self-overlap method (SO) [lE El, 
has already been successfully used in the study of the 
stability and critical points of random Boolean networks, 
a system that is multi-component albeit deterministic. 
In this work we will show that the method can also be 
successfully applied to a stochastic system such as a 
spin network. We will obtain analytical and numerical 
results on the 2D Glauber-Ising model that exactly 
match those yielded by DS. However, contrary to DS, 
SO proceeds by handling only one replica and analyzing 
its own evolution in time -using basically one-point 
correlation functions at subsequent time steps for that 
task. The computational costs are thus lower in SO 
than in DS. As we will see the analytical calculations 
also become much simpler while yielding the same 
results. Furthermore, and what is more important, the 
SO method is free from the ambiguities that afflicted 
DS due to its use of two replicas. This comes as a 
direct consequence of the already mentioned fact that 
on a single replica it does not matter whether one uses 
Glauber or HB dynamics since they posses the same 
one-point correlation functions. 

We will follow the development applied by Vojta in 
, comparing in each case the results obtained using DS 
and our results (using SO). The paper is organized as it 
follows: in section II we quickly introduce both the 2D 
Glauber-Ising model and SO. We then apply the method 
to the 2D Glauber-Ising model in section III, obtaining 



a system of equations (master equation) that describe 
the dynamical evolution of the system. We discuss then 
how to apply a mean-field approximation to the system, 
and compare it with the methodology used by Vojta 
12]. In section IV we obtain an analytical expression 
for the magnetization of the system in both ferromag- 
netic/paramagnetic phases similar to that obtained by 
Vojta [2] . Numerical results are provided at this point in 
order to validate the mean field approximation assumed 
in the analytical development. Finally, in section V we 
provide some analytical and numerical results on the sta- 
bility of the model, showing that the system is chaotic 
(disordered) in the paramagnetic phase. Conclusions are 
presented in section VI. 



II. ISING MODEL, THE DAMAGE SPREADING 
VS. SELF-OVERLAP METHOD 

A. Glauber Ising model 

We will work with a kinetic Ising model, a lattice of 
N spins. Si e {+1,-1}, that follows Glauber dynamics. 
That is, at every time step a lattice site i is chosen at 
random. If the spin value of site i at time t is given by 
Si(t), at time t + I it will be given by: 



Si{t + 1) = sgn 



(1) 

where £,i{t) is a random number such that S,i{t) G [0, 1). 
The transition probability (/>(ft.i) is given by the usual 
Glauber expression: 



(2) 



where T denotes the temperature and hi{t) is the local 
field seen by spin i at time t: 



j—n.n 



JijSj{t) + ho. 



(3) 



In this expression Hq represents an external magnetic 
field, and the sum in the interaction term applies only 
to the nearest neighbors (three for example in an hexag- 
onal lattice) . Without loss of generality from now on we 
will take /iq = and Jij = 1. 



B. Damage spreading and self-overlap 

As stated above we will use the SO method to study 
the dynamics of the system. This procedure was intro- 
duced by Luque and Ferrera [3] and its underlying phi- 
losophy is similar to that of the DS method used by Vo- 
jta to study the thermodynamics of phase transitions in 
spin systems. The main difference between the two pro- 
cedures lies in that, while damage spreading uses two 
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copies of a system with slightly different initial condi- 
tions (the damage) and computes the evolution of these 
differences, the self-overlap method uses the difference 
between successive temporal states of a single system as 
the system evolves towards equilibrium. For instance, in 
DS the damage D{t) at time t is defined as: 



— T 

2N ^ 



,(1) 



,(2) 



(4) 



and measures the (averaged) Hamming distance be- 
tween the states of the two replicas at that time (i.e., 
the proportion of sites for which the spin state differs be- 
tween the system (1) and the damaged replica (2)). In 
SO however the self-overlap a{t) at time t is defined as 
one minus the averaged Hamming distance between the 
states of a spin site at time t — \ and at time t: 



a(t) = 1 - — V 

i=l 



S,{t) - S,{t - 1) 



(5) 



In order to describe the time evolution of the system 
it is useful to define the "up state self-overlap" a++(t) 
at time t as the average number of spin sites that had 
Si — +1 both at time t — 1 and at time t. We also 

define a (t), (t) and a ^(t) in a completely similar 

fashion. By normalization we must then have: 



a++(t) + a__(t) + a+_(t) + a_+(t) = 1. 



(6) 



Since the sites that remain in the same state at times 
t — 1 and t drop from the sum in the definition ^ we 
also must have 

a{t) ^ 1 - a+_(t) - a_+(t) = a++{t) + a__(t). (7) 

Once the equilibrium as been reached the relation = 

a \. must be satisfied, where we have dropped the time 

dependence to indicate equilibrium values. Then trivially 
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(8) 



At this point it is interesting to note that the self-overlap 
functions can be understood in terms of autocorrelation 
functions, more precisely, two-time autocorrelation func- 
tions. For instance, in equation ([5]), one can rewrite 
\si{t) — Si{t — 1)1 as (1 — Si{t)si{t — l))/2, which is a 
shifted autoresponse function measured at subsequent 
time steps. In a similar way, the rest of self-overlap func- 
tions can be written as linear combinations of the basis 
of autoresponse functions. 

Autocorrelation functions have been widely used as effi- 
cient tools in order to measure spatial or temporal corre- 
lations in physical and biological systems (repeated pat- 
terns, relaxation, frustration, etc). Their applications 
ran ge f rom investigations in transport properties of flu- 
ids [2l| or the analysis of climatological models [23] to 



studies of decoherence in quantum systems ^3|, to cite 
but a few. Autocorrelation functions are the center of 
interest in theoretical studies of the relaxation of non- 
equilibrium systems. In this sense, much work has been 
recently done in order to characterize dynamical scal- 
ing and other invariant behavior in the ageing regimes 
of Ising-like systems [13, IM HI, 113, HI ■ In our case it 
would be fair to say that the self-overlap functions are re- 
ally measurements of autocorrelations under a different 
garment. To dwell on a deeper review of the existing lit- 
erature on autoresponse functions would go beyond the 
scope of this paper however. We would like to emphasize 
nonetheless that what is new here is: i) the fact that this 
particular combination of self-correlation functions mea- 
sured at subsequent time steps manages to capture the 
essence of the (same-site) temporal correlations in sys- 
tems that undergo order/disorder phase transitions, and 
ii) this is then combined with a philosophy inspired by 
DS, namely: an evolution equation towards the equilib- 
rium state for the correlations, and a mean field approxi- 
mation directly extracted from DS in order to be able to 
solve this equation. Once the evolution equation and the 
mean field approximation are in place the self-overlaps 
will allow us to study the stability of the different states 
accessible to the system, and hence the phase transition 
itself. 



III. MASTER EQUATION, TRANSITION 
PROBABILITIES, AND MEAN FIELD 

A. Master equation 

Generally speaking, the self-overlap method would pro- 
ceed by solving some master evolution equation for the 
a's in order to obtain their equilibrium values, much in 
the vein of the damage spread method. We begin by 
defining the probability of finding a spin site in the -|- 
(-) state at time t, P+{t) {P-{t)) 



P±it) 



n±{t) 
N 



(9) 



where n±(t) is the number of sites with spin up (down) 
at time Obviously P4.(t)-|-P_(i) = 1. By the definition 
of a++(i), a \-{t) it follows that 

P,it)=a^,it) + a,,it)^l±^^±±^^^^ (10) 

and analogously for the down states 

P_(t)^a,_(t)+«__(t)^ ^ + "-(^)-"^+(^) . (11) 

As noted above in the limit t og the a's ought to 
reach their equilibrium values and one can drop the t 
dependence. 

Of particular interest to us will be the transition prob- 
abilities from one state to another, i.e., the elements of 
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the transition matrix of our Markov process. Let W++{t) 
{W^^(t)) be the average probabihty of changing from the 
+ (— ) state at time t to the + (— ) state at time t + 1, 
where the precise meaning of this average wih be made 
clear shortly. In an mean field approximation we will 
then have 



a++it)^W++it-l)P+it-l), 



(12) 



and analogously with ,W |-. Note that these W^s 

will then be the elements of an average Markov matrix 
for the evolution of the system. Combining (jlO|l .(fTT |) and 
together it is easy to arrive to a couple of mean field 
evolution equations for 0++ (t) and a (t) , namely 



M+ + -} = m 



,1/T 



2cosh(l/T)' 



-l/T 



,{ } = <^(-3) = 



2 cosh (-l/T)' 
2cosh(-3/r)' 



(15) 



Note that the calculations are much simpler than those 
needed in DS 0. The probability associated to each 
configuration will be 



p3 



P(01) = 3P2(1-P+ 



-a++(t) = -a++{t)W+-{t)+a^+{t)W++{t), 



dt 



a__(i) = -a^^(t)W-+{t) + a+-{t)W—{t).{13) 



These two equations are of course nothing but the 
reaction-diffusion equations for the a's that common 
sense would have dictated us to begin with. We now 
proceed to evaluate a mean field approximation for the 
VF's so that we may solve 



P(02)=3P+(1 



Pi 



(16) 



where to simplify the notation we have dropped the time 
dependence, although in this case one must be aware that 
we are not dealing with equilibrium values (this will be 
the case for the next several equations). Using equations 
pop and pi[) . we can now write after some trivial ma- 
nipulations 



B. Mean- field approximation 

To begin with, note that in a system that follows 
Glauber dynamics the transition probability at site i for 
a given local field hi is given by ([2]) above. This means 
that 



W.+ih,) = 0(/z,), 



W+^{h,) = l~(t>ih,), 



(14) 



That is, as is well known for a given local field hi the 
probability that the spin at site i will be in the -I- state 
at time 1 is always 4>{hi), whereas the probability that 
its state be — will be 1 — <t>{hi), regardless of the initial 
state of the site. Thus, finding average values for the W^'s 
is equivalent to finding an average 4>{hi), (j). 

The mean field approximation that we will use closely 
follows the spirit of the effective-field approximation used 
by Vojta [2,]. This consists basically in averaging over all 
the possible configurations that can surround a given site, 
where in the average each configuration is weighted by 
its probability of taking place. Thus, with three nearest 
neighbors per site, the transition probabilities can take 
the values (remember that we are taking Jij = 1) 



0o{+ + +} = 0(3) 



c3/T 



2cosh(3/r)' 



k=0 



1 3 , 

2 + 8(«++-«--) 




(17) 



Using the relations between the a's and applying the 
mean field to the right hand side of the differential equa- 
tions p3p we can rewrite them as 



dt 



a++ 



-a++(l 



1 — a++ — a 



^ (18) 



dt 



-a = — a_ 



1 — — a_ 



(1 -(/.), (19) 



which by P7)) is now a system of equations depending 

only on a-|-+ and a Note that it is easy to generalize 

the mean field approximation to the case of n nearest 
neighbors (that is, for a given topology): 



^ \k 

fc=0 



pn — kpk 



1 



1 



exp 



(20) 



This would be much harder to do using DS, if at all pos- 
sible. 
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IV. THERMODYNAMICS: MAGNETIZATION of configurations that need to be taken into account. 



At this point we are going to link the self-overlaps 
to the average magnetization per spin, m. With 
P+ (t) , P_ {t) as defined above 



P±{t) 



n±it) 
N ' 



(21) 



we must then obviously have for the average magnetiza- 
tion m 



m(i) =P+(i)-P_(t) 
or, since P+(t) +P_(t) = 1, 
1 + m(t) 



(22) 



„ , , 1 — m(t) 
P-it) = (23) 



2 ' ' ' 2 

By the definition of a++ (t) , a ^(t) it follows then 

a_+(i) + a++(t) = P+{t) = i±|^, (24) 

and analogously with a__(i),a-| (t) and P-{t). 

Since 

TO(i) = a++(i) - a__(i), (25) 
the system of equations (|17I18I19P can be rewritten as 



d TO f ^ 3 
— TO = — <^ -1 + - 
dt 2 1 4 



tanh 



tanh ( — 
T 



tanh ( — 



T 



3 tanh 



T 



(26) 



Within the limits of our approximation this equation de- 
scribes the evolution towards equilibrium of the magne- 
tization TO for the case of n = 3 nearest neighbors. Set- 
ting dm/dt = one can obtain an expression for the 
temperature dependence of its equilibrium value ra{T), 
and from it one can extract the transition temperature 
for the ferro-paramagnetic transition — this was the ap- 
proach originally followed by Vojta \2\. 
Equation yields a critical temperature Tc « 2.104 
above which the magnetization is zero. When T < Tc, 
we have: 



m — ±\ 



1 + f [tanh (f) + tanh (^)] 
I tanh (i) - i tanh (^) 



(27) 



Both results completely coincides with those in . Note 
however that the calculations involved here have been 
considerably simpler — again basically due to the fact 
that in SO we only consider one replica of the system, 
which results in a considerable reduction in the number 



In our Monte Carlo simulations, the procedure to mea- 
sure the (equilibrium) self-overlap goes as follows: let 
us suppose that we generate a random initial condition 
for the N spin lattice. Then we let it evolve towards 
equilibrium by applying the Glauber dynamics with 4 
neighbors (square lattice). Once equilibrium has been 
reached we compute the states of the system for a suf- 
ficiently large number of time steps. We have used in 
all cases 10, 000 x N time steps for a square lattice of 
N — 100 X 100 spins (that is, defining a system time step 
t sls N steps of the simulation, we use t = 10000 system 
time steps). If we then count the number of times that 
a spin site is in the "up" state, both at time t and 
t — 1 and average over all sites and time steps, this will 
give us the equilibrium value of the up state self-overlap 
a++. Repeating this procedure with the down state, — , 
will then obviously give us the "down state self-overlap" , 

a , and so on. Each value of the simulation is averaged 

over 100 realizations. 

In figure ([1]) we plot the average equilibrium magnetiza- 
tion vs. temperature in order to visualize how our mean- 
field approximation performs — we note here that we are 
basically interested in the thermodynamic limit of infi- 
nite lattice size and that we are removing the inherent 
degeneracy of the system by plotting only positive mag- 
netization. First, note that our Monte Carlo simulations 
in a square lattice (squares) are in fair agreement with 
the Onsager (infinite size) solution -dashed line- except 
in the proximity of the phase transition, where finite size 
effects are relevant and difficult to suppress. Compar- 
ing then the Monte Carlo simulations and the mean-field 
solution with n — A neighbors we can see that qualita- 
tively speaking they provide the same results, with the 
mean field typically overestimating the critical temper- 
ature. We stress here however that the purpose of this 
paper was not so much to present a mean field technique 
able to reproduce the exact results, but rather to intro- 
duce a new technique able to exactly reproduce previ- 
ously known mean field results while at a much lower 
cost. For illustrative purposes and to allow comparison 
with the results obtained by Vojta we also show in figure 
[1] the mean-field result for n = 3 neighbors (hexagonal 
lattice), which underestimates the n = 4 critical temper- 
ature Tc- 

In figure ([2]) we plot the equilibrium values , a* 

vs. temperature, following the same methodology of 
figure ll]): we compare our Monte Carlo simulations 
(circles) with the numerical resolution of the mean field 
equations (note that again, the mean field with n — 3 
underestimates the quantitative behavior and the one 
with n = 4 overestimates it). As we can see in the 

figure, the self-overlap a = a++ + a acts as an order 

parameter. 

We also note that more work remains to be done in 
order to make an in-depth comparison between DS and 



6 




FIG. 1: Magnetization of the system versus temperature in 
the case of: (squares) Monte Carlo simulation of a 100 x 
100 spin square lattice (the solid line here is just a guide for 
the eye), with 10,000 system steps and averaged over 100 
realizations ; (dashed-dot line) mean field approximation for 
n — 3 neighbors; (solid line) mean field approximation for 
n = 4 neighbors; (dashed line) Onsager solution. Note that 
the mean field approximation recovers the expected behavior, 
that is, null magnetization above Tc, non null magnetization 
below Tc, which tends to a constant maximum value at T = 0. 
The difference lies on the quantitative value of Tc in each case, 
overestimated by the mean field in the case of n = 4 neighbors. 

SO. For instance, one may evaluate the critical expo- 
nents of the self-overlap order parameter a and compare 
the results with the DS approach _9] , which would be in- 
teresting. This however goes somewhat beyond the scope 
of this paper. 

V. STABILITY 

In equation ([26]) , the stabihty of the fixed point m* = 
(paramagnetic phase) is related to the sign of the eigen- 
value: 

A(T) = l/2(-l + 3/4[tanh(l/T) + tanh(3/r)]). (28) 

Note that (|26p falls into the normal form of a pitchfork 
bifurcation at T = Tc where the fixed point is not hy- 
perbolic and the Hartman-Grobman theorem [2^ does 
not apply. For T > Tc, m* = is stable, and below it, 
it becomes unstable. The fact that we have a pitchfork 
bifurcation at Tc implies that in the ferromagnetic phase 
(i.e., below Tc) two other stable fixed points must appear. 
They are indeed ±m* , where m* is now given by (|27p . 
Taking into account the relation between m and a, with 
a little algebra we arrive at 

a* =m*(l-^*) + 0*. (29) 

Hence, the fixed point m* — leads to 0* = 1/2 (ac- 
cording to (fT7|) ) and a* — 1/2, which are thus stable at 



FIG. 2: Stationary values of a++ and a in the case of: 

mean-field approximation with n — 3 first neighbors (dashed 
line), mean-field approximation with n — 4 first neighbors 
(solid line) and Monte Carlo simulation of a 100 x 100 spin 
square lattice (here the solid line is just a guide for the eye), 
with 10, 000 system steps and averaged over 100 realizations 
(circles). Note that at Tc a pitchfork bifurcation takes place 
in the three cases. The bifurcation value is underestimated 
by the mean field approximation in the case of n = 3 neigh- 
bors and overestimated in the case of n = 4 neighbors. Below 

the critical temperature a++ — > 1 while a ^0 for a system 

that chooses the m — +1 vacuum, whereas the opposite is true 
if the system goes to m = — 1. Above Tc the system tends to 
(a++,a ) = (l/4, 1/4). Note that although the critical tem- 
perature is only predicted qualitatively, the stationary values 
for (a++,a__) yielded by our simple model exactly match the 
Onsager predictions. 

T > Tc- Note that a = 1/2 is the minimal self-overlap 
that the system can show. 

We can define at this point a Hamming-like distance be- 
tween successive temporal states (a self-distance), as: 

d{t) = l-a{t). (30) 

The fixed point a* = 1/2 implies that we must have a 
fixed point for d at d* = 1/2 which, since it is taking 
place at the minimal self-overlap, is equivalent to the 
maximal self-distance of the system (total disorder). 
Following Wolf's method as in the case of random 
Boolean networks [l^, this self-distance would enable 
us to determine a Lyapunov exponent of the system. 
However, one can simply apply the Hartman-Grobman 
theorem directly [1^. Near the fixed points the self- 
distance of our system can be expressed in terms of 
d{t) ^ exp(At), where A is given by (|28p . This eigenvalue 
can also be understood as a Lyapunov exponent. Note 
that nevertheless it would not be a standard Lyapunov 
exponent: when d{t) tends to its fixed point, the system 
is actually tending to the maximal disorder, thus A < 
means chaos. 

Summing up, in the paramagnetic phase, m* = is 
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FIG. 3: Values of the temperature dependent eigenvalue (|28p 

of J: when it is negative, (a++,a ) = (l/4, 1/4) is stable, 

thus the self-distance d = 1/2 is the attractor of the system 
(chaotic phase). When the eigenvalue is positive, the value 
(1/4, 1/4) is a saddle point and thus an unstable fixed point. 
At Tc ~ 2.104 the eigenvalue is null, thus the fixed point is 
not hyperbolic -a bifurcation takes place-. 

stable, thus d* — 1/2 is stable too: the system tends 
exponentially to the maximal disorder and the phase is 
chaotic. 

Figure ([2]) is a plot of equation ([^5]) . Note that when 
T > Tc (paramagnetic phase) an increase of the temper- 
ature leads to an increase of chaos, with the self-distance 
of the system tending faster to the attractor d* — 1/2. 

In the ferromagnetic phase however the stable sta- 
tionary value of d is 

d*(T) = l-m*(T)(l-0*)-0*, (31) 

with m* given by (j27p and the fixed point value of the 
mean field. The self-distance tends to zero for low T, and 
thus the system is in a frozen state (order) . When we in- 
crease the temperature the self-distance also increases up 
to the maximum value d — 1/2, which is reached at Tc 
(no correlation) . These results agree with those found in 
the paramagnetic phase. We can conclude therefore that 
our approach correctly reproduces an ordered behavior 
in the ferromagnetic phase and disordered (chaotic) be- 
havior in the paramagnetic phase. In Appendix (A) we 
perform a more detailed analysis of the stability of the 
system that confirms this conclusion. 

VI. CONCLUSION 

In this paper we have introduced the self-overlap 
method by using it to study both analytically and nu- 
merically the 2D Ising model. Since the properties of this 
model are obviously well known our main concern was to 



show that SO is an unambiguous method (with respect 
to changes in the algorithm implementation) that cor- 
rectly reproduces the standard results while being very 
advantageous from both the numerical and the analytical 
point of view. The SO method could thus constitute a 
rather simple and eSicient method of stability analysis in 
this kind of multicomponent systems (Ising-like models, 
spin glasses, CA, Kauffman networks, etc). Many other 
physically relevant quantities in these systems (measures 
of complexity, information theory measures such as the 
mutual information, and so on) can be studied and mea- 
sured by applying SO, something that we think deserves 
further investigation. Wherever damage spreading was 
supposed to have been useful and the equilibrium state 
of the system is ergodic, we think that self-overlap ought 
to work too and do so in a non ambiguous manner. More- 
over, it should also be more efficient numerically speak- 
ing, and simpler from the analytical viewpoint. 
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APPENDIX A: DETAILED ANALYSIS OF THE 
STABILITY 

We undertake here a deeper study on the stability of 
the system. For this task we go back to the evolution 
equations [TO]) , which constitute a nonlinear differen- 
tial system. The fixed points of this system are obtained 
from equating (|18ll9p to zero (reducing the differential 
system to a linear system). This yields a total of three 
fixed points, namely: 

(a++(T)*, a__(T)*) = (^|(1 + m*), m*(| - 1) + , 
(a++(r)*, a__(T)*) = (m*i^ _ l) + ||(i + „,*)^^i) 

when T < Tc (where m* is given by (|27p ). and 
(a++(T)*,a__(T)*) = (1/4,1/4) VT (this solution is ob- 
viously related to the fixed point m* — 0). 
We can write (f> as 

0= ^ + ^^(^)(a++-a-) + ii?(^)(a++-a__)^ (A2) 
where 

A{T) = ^[tanh(3/r) + tanh(l/T)], (A3) 

and 

B{T) = [tanh(3/T) - 3tanh(l/r)]. (A4) 
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Let's start with the stabihty analysis of the fixed point 
(a++,al„) = (1/4,1/4). This solution is independent 
of T and for T > Tc is the only fixed point (note that 
in this case </) takes the value 1 /2 independently of the 
number n of neighbors as it can be proved after some 
trivial algebra) . Computing the jacobian J at this fixed 
point, we come to: 

■J l(i/4,i/4)-4 A{T) - 1 -A(T) -3 )' 

with eigenvalues Ai = -1, and A2 = 1/2(A(T) - 1). We 
will distinguish then three situations: when A(T) < 1, 
(1/4,1/4) is an hyperbolic (indeed stable) fixed point 
(which is obviously related to the fact that m* — is 
stable when T > Tc). When A{T) > 1 the fixed point is 
again hyperbolic, but now it is unstable (a saddle point). 
In these two situations we can apply the developed 
formalism, due to the Hartman-Grobman theorem [29|. 
Hence, A{T) < 1 ^ T > 2/ln(22/3 + 1) w 2.104 (and 
viceversa for A{T) > 1). 

We thus get that when T > Tc (that is, in the 
paramagnetic phase), the stationary solution (1/4,1/4) 
is stable. In the ferromagnetic phase however (T < Tc) 
this fixed point becomes unstable. 

At this point we can introduce the self-distance 
defined in The stability of the (1/4, 1/4) solution 

directly implies that d will have a stable value of 1/2 in 
the paramagnetic phase, whilst this value will become 
unstable in the ferromagnetic phase. Since in the 
paramagnetic phase (1/4,1/4) is the only fixed point 
the self-distance necessarily goes to the attractor (stable 
fixed point) d* = 1/2, indeed exponentially due to 
the Hartman-Grobman theorem, and the phase is thus 
chaotic. However in the ferromagnetic phase (1/4, 1/4) is 
unstable: orbits with initial conditions arbitrarily close 
from this fixed point will separate from it exponentially, 
correlations will take place and the phase will become 
ordered. 

When A{T) = 1, applying Peixoto's theorem [2^ . 



we can conclude that (1/4,1/4) is a bifurcation point 
(lack of structural stability), that is, Tc constitutes a 
bifurcation value. What kind of bifurcation is taking 
place?. It is easy to see that the linearized system has a 

symmetry of the type a++ — a Using this symmetry, 

the system of equations (|18ll9p can be transformed into 
(PSl) . This equation falls into the normal form of a codi- 
mension one bifurcation, a pitchfork bifurcation (indeed, 
subcritical). This means that two branches of equilibria 
appear for T < Tc associated with values of m 7^ 0, 
either positive (positive branch) or negative. Undoing 
the change of variables we get that below Tc we must 
have, for a given T, two extra stationary points -other 
than (1/4,1/4)- of the shape [(a, 6), (6, a)]. These fixed 
points correspond obviously to (jAip . Moreover, since as 
the Poincare index is a topological invariant these two 
new fixed points are both stable in the ferromagnetic 
phase (in the paramagnetic phase the global index is -f 1 
because the fixed point (1/4,1/4) is a sink, whereas in 
the ferromagnetic phase (1/4, 1/4) is a saddle point with 
index —1, so the other two fixed points must have index 
+1). Depending on the initial conditions, the system will 
evolve to a fixed point of the shape (a, 6) or to (&, a). In 
other words, the Ising model will give us either positive 
or negative magnetization in the ferromagnetic phase, 
depending on the initial condition. If the system starts 
at T > Tc, where the magnetization is zero, and we 
lower its temperature below the critical one, fluctuations 
will take the system either to the upper or to the lower 
branch indistinctively. 

In figure ^ we plot together the stationary values 
{a'^_^,a*__) of the differential system (|18ll9p for both 
n = 3 and n — A nearest neighbors and the results 
from our Monte-Carlo simulation (again, a square lat- 
tice of 100 X 100 spins, where we ran 10, 000 system steps 
after reaching equilibrium, and averaging over 100 real- 
izations). We can see that the results are qualitatively 
similar, that is, the stationary value (1/4, 1/4) is stable 
above the Curie temperature and unstable below it. As 
expected, at Tc a pitchfork bifurcation takes place and 
when T <Tc the system has two stable fixed points, i.e. 
(a, b) and (6, a) for each T. 
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